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ABSTRACT 

We  present  a  continuation  technique  for  branches  of  periodic  solutions 
which  can  be  applied  to  autonomous  systems  of  nonlinear  differential 
equations.  Coupling  the  technique  with  Hopf  bifurcation  expansions  enables 
one  to  compute  entire  periodic  solution  branches  including  those  with  turning 
points  and  unstable  solutions.  We  apply  these  methods  to  a  continuously 
stirred  tank  reactor  with  consecutive  A  +  B  +  C  reactions.  Our  computation 
reveal  dynamic  phenomena  not  seen  in  previous  reactor  studies.  The  results 
include  response  diagrams  exhibiting  stable  and  unstable  periodic  branches 
that  contain  multiple  turning  points.  The  discovery  of  these  points  indicate 
the  reactor  may  jump  from  a  steady  state  to  a  periodic  orbit  or  may  jump  from 
one  periodic  orbit  to  anther. 
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SIGNIFICANCE  AND  EXPLANATION 


Numerical  methods  are  developed  for  tracing  branches  of  periodic 
solutions  exhibited  by  mathematical  methods  consisting  of  autonomous, 
nonlinear  differential  equations.  These  continuation  methods  allow  one  to 
compute  the  periodic  orbits  which  emanate  from  steady  state  solutions  at  the 
Hopf  bifurcation  points.  It  should  be  emphasized  that  the  techniques  permit 
us  to  compute  entire  periodic  branches  including  those  with  turning  points, 
unstable  solutions  and  infinite  periodicity. 

We  apply  the  methods  to  a  mathematical  model  of  a  stirred  tank  reactor 
with  an  A  ♦  B  ♦  C  reaction  mechanism.  The  reactor  exhibits  broad  regions  of 
multiple  steady  states  and  periodic  solutions,  and  we  discuss  the  effects  of 
the  intricate  solution  interactions  on  reactor  operation. 
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NUMERICAL  COMPUTATION  OF  PERIODIC  SOLUTION  BRANCHES  AND  OSCILLATORY 
DYNAMICS  OF  THE  STIRRED  TANK  REACTOR  WITH  A  ♦  B  ♦  C  REACTIONS 

Eusebius  J.  Doedel*  and  Robert  F.  Heinemann 

INTRODUCTION 

The  bifurcation  analysis  of  time-dependent  differential  equations  has  found 
application  in  several  fields  of  research  including  chemical  reaction  engineering  [1-8], 
mathematical  biology  and  biochemistry  [9-16],  and  combustion  [17-24],  The  first  step  in 
analyzing  the  mathematical  models  from  these  fields  normally  consists  of  determining  all 
the  possible  steady  state  solutions.  In  general,  such  computations  require  a  continuation 
procedure  for  steady  state  solution  branches  which  allows  for  computing  past  turning  points 
without  difficulty  and  for  switching  from  one  branch  to  another  at  points  where  the 
branches  intersect.  Techniques  developed  by  Keller  [25]  fulfill  these  requirements  and  are 
used  in  studies  of  multiple  solutions  exhibited  by  tubular  reactors  [26],  catalyst  pellets 
[27],  premixed  flames  [20-22]  and  flow  between  rotating  cylinders  and  disks  [28,  29], 

The  second  step  in  the  model  analysis  is  often  the  determination  of  the  periodic 
solutions.  When  a  branch  of  such  solutions  crosses  a  steady  state  branch,  the  point  of 
intersection  is  called  a  Hopf  bifurcation  point.  In  reaction  engineering,  Hopf  theory  has 
been  used  to  determine  all  possible  periodic  orbits  exhibited  by  an  exothermic  A  ♦  B 
reaction  in  a  continuously  stirred  tank  reactor  (CSTR)  [6,  7].  This  local  theory  provides 
perturbations 1  expansions  at  the  Hopf  point  which  determine  the  frequency,  direction, 
stability  and  profile  of  the  bifurcating  solution  without  time  integration.  The 
application  of  the  expansions  is  particularly  useful  in  examining  reactor  dynamics  and  jump 
phenomena  which  accompany  the  stability  exchange  at  the  Hopf  points.  Recently  these 
techniques  were  extended  to  distributed  parameter  systems  and  applied  to  the  tubular 
reactor  [8]. 
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Further  insight  into  the  dynamics  of  chemical  reactors  is  achieved  by  tracing  the 
periodic  solution  branches  away  from  the  bifurcation  points.  These  computations  are 
usually  accomplished  by  applying  initial  value  techniques  to  the  model  equations  and 
performing  transient  calculations  until  the  solutions  converge  into  the  periodic  orbit. 

This  strategy  is  capable  of  determining  stable  periodic  solutions  as  well  as  unstable 
solutions  for  the  special  two  dimensional  case  (e.g.  the  CSTR  with  an  A  ♦  B  reaction) 
where  one  can  simply  replace  t  by  -t  in  the  model  equations.  In  most  other  cases, 
initial  value  techniques  will  be  difficult  or  impossible  to  use,  especially  when  turning 
points  occur  along  the  periodic  branches  (30). 

In  this  paper  we  present  a  simple  continuation  method  for  periodic  solutions  which 
allows  the  computations  to  proceed  past  such  turning  points  as  well  as  along  unstable 
branches.  We  transform  the  initial  value  problem  into  a  boundary  value  problem  on  the 
fixed  interval  (0,2* ]  by  treating  the  period  of  the  solution  as  an  unknown  and  imposing 
periodicity  as  a  side  condition.  Two  equations  are  added  to  regulate  the  step  sire  along 
the  solution  branch  and  to  avoid  translation  of  the  periodic  solution.  The  resulting 
system  of  coupled,  nonlinear  equations  is  efficiently  solved  by  Newton's  method.  A  similar 
treatment  for  tracing  branches  of  periodic  solutions  has  bean  used  by  Rinzel  and  Miller 
[10]  and  has  received  extensive  treatment  in  a  recent  thesis  [31].  The  particular 
formulation  of  the  numerical  scheme  that  we  use  however  has  great  practical  advantages 
[30].  Other  numerical  techniques  for  Hopf  bifurcation  are  reported  by  Weber  [32]  and 
Langford  [33]. 

We  apply  our  numerical  techniques  to  an  important  three  dimensional  problem  - 
consecutive  A  *  B  ♦  C  exothermic  reactions  in  a  CSTR.  In  an  earlier  paper,  Hlavacek, 
Kubicek,  and  Visnak  [34]  give  a  rather  complete  classification  of  the  steady  state 
multiplicity  for  this  system  and  also  present  isolated  numerical  calculations  illustrating 
the  existence  of  oscillatory  solutions.  Responding  to  the  need  pointed  out  by  Hlavacek 
et  al.  for  work  in  multi-dimensional  systems  exhibiting  oscillations,  Cohen  and  Keener  [35] 
apply  multi-time  scale  perturbation  techniques  and  illustrate  the  utility  of  this  method  by 
deriving  formulas  for  the  bifurcating  solutions  for  a  few  combinations  of  model 


parameters.  However,  the  most  thorough  classification  of  dynamics  for  this  reactor  to  date 


is  given  by  Halbe  and  Poore  [36] .  They  apply  the  Hopf  theory  to  illustrate  the  increased 
diversity  of  the  consecutive  reaction  case  compared  to  the  A  ♦  B  example  and  show  that 
the  A  *  B  ♦  C  CSTR  may  exhibit  up  to  five  steady  state  states  and  up  to  four  Hopf 
bifurcation  points. 

None  of  the  above  studies  traces  the  periodic  branches  away  from  the  Hopf  point.  We 
present  in  a  later  section  entire  solution  branches,  containing  stable  and  unstable  orbits, 
which  are  computed  using  our  numerical  techniques.  Although  we  do  not  intend  to  provide  a 
complete  cataloguing  of  the  dynamic  capabilities  of  this  reactor,  our  computations  reveal 
phenomena  not  reported  in  earlier  CSTR  studies.  Of  particular  Interest  are  branches  which 
contain  multiple  turning  points  creating  the  possibility  of  jump  transitions  between 
periodic  states. 


MASS  AMD  ENERGY  BALANCES 

we  consider  consecutive  exothermic,  first  order  reactions  in  an  ideal  stirred  tank 
continuously  fed  by  component  A  and  cooled  by  a  medium  with  constant  temperature  Tc. 
The  mass  balances  for  A  and  B  and  the  energy  balance  are  written  below  [34] 

(1) 

(2) 

PC  V  7T  -  PC  v(Tn  “  T)  "  US  (T  -  T  )  +■  (-4b,  )VC,A,exp(-E,/RT)  (3) 

pat  p  o  a  c  i  A  i  1 

♦  <-dH2)VCBA2exp(-E2/RT) 


v  dT  '  V(CA0  -  CA>  -  VCAA1exp(-E1/RT) 
dCB 

V  dt“  "  _VCB  +  VCAA1exp(-E1/RT)  -  VCBA2exp(-K2/Rt) 


These  equations  can  be  written  in  dimensionless  form 


dj- 

dT 


1  -  y  -  Dy  exp(8/1  +  e8) 


(4) 


~  «  -z  +  Dy  exp(8/|  +  c8)  -  DCz  exp(1  *9e9)  (5) 

,0  _  g 

—  »  -8  -  8(8  -  8^)  +  DBy  exp(8/l  +  e8)  +  DBao  exp( ^  ■— (6) 

Here  y  is  dimensionless  concentration  of  X,  z  is  dimensionless  concentration  of  B, 

8  is  dimensionless  temperature,  and  T  is  dimensionless  time.  D  is  Damkohler  number, 
e  is  dimensionless  activation  number,  0  is  selectivity  ratio,  K  is  the  ratio  of 
activation  energies,  8  is  the  dimensionless  heat  transfer  coefficient,  B  is  the 
adiabatic  temperature  rise  and  a  is  ratio  of  heats  of  reaction.  It  is  convenient  to 
analyze  the  limiting  case  of  small  e  and  8^  -  o. 

NUMERICAL  METHODS 

The  numerical  methods  presented  in  this  section  are  capable  of  treating  nonlinear 
autonomous  systems  of  ordinary  differential  equations: 

-  f(u(t)»X),  t  >  0,  u,f  B  Rn  ,  (7) 

where  X  is  a  free  parameter.  A  more  detailed  presentation  of  these  numerical  procedures 
and  computer  codes  can  be  found  in  130] . 

Steady  State  Computations 

The  computation  of  the  steady  state  solution  branches  is  an  algebraic  problem 
consisting  of  the  bifurcation  analysis  of  the  equation 

f(u,X)  »  0  .  (8) 

This  is  accomplished  numerically  by  applying  the  continuation  and  switching  techniques  due 
to  Keller  (25).  This  method  inflates  the  original  problem  by  imposing  an  arc- length 
normalization,  N,  on  the  solution  branch  so  that 


f (u(s) iX(s) ) 
N(u( s) ,X( s) >  s) 


(9) 


where  the  new  unknown  s  (a  "pseudo-arc-length"'  replaces  X  as  the  continuation 
parameter  in  an  Euler-Newton  scheme.  The  main  advantage  of  using  the  inflated  system  is 
the  capability  of  computing  past  turning  points  on  the  solution  branch,  and  this  idea 
serves  as  a  basis  for  the  continuation  techniques  used  to  trace  the  periodic  solutions. 
Periodic  Solution  Branches 

In  computing  entire  branches  of  periodic  solutions  for  autonomous  systems  both  the 
solution  u(t)  and  the  solution  period  p  continuously  change  as  X  varies.  To  fix  the 
period  we  scale  the  time  variable  t  and  transform  (7)  into 


XT  -  fr  f(u(t)»X)  t  e  [0,2s] 


(10) 


Now  the  unknown  period  p(X)  appears  explicitly  and  we  seek  solutions  satisfying 

u(0)  -  u(2»)  (11) 

Thus  the  initial  value  problem  has  been  replaced  by  a  boundary  value  problem  with  non- 
separated  boundary  conditions  which  allows  the  computation  of  unstable  orbits.  Suppose  we 
are  given  a  periodic  solution  (v(t) ,P0,Xg).  Our  objective  is  to  find  a  neighboring 
solution  (u(t),P,X)  at  some  distance  As  from  the  given  solution.  To  this  end,  we  solve 


(12) 


where  C  and  B  represent  (10)  and  (11).  The  unknown  period  P  is  included  in  the 
normalisation 


-5- 


0 


(13) 


N(U,P,X)  =  62lu  -  V|2  +  02,x  _  XQ|  +  83(P  -  P0)2  -  is2  - 


where  6 ^ *s  are  fixed  weights  and  ***2  **2  norm* 

!  A  represents  an  "anchor"  equation  and  is  included  to  insure  unique  periodic  solutions 

j 

by  eliminating  free  translation  in  time.  Although  there  are  many  possible  choices  for  an 
appropriate  anchor  equation,  the  following  is  found  to  be  extremely  convenient  for  problems 
with  rapidly  varying  solution  components 

i 

2* 

Mu)  S  /  (u(t)  -  v(t))  f(v(t)ji  )dt  »  0  (14) 

0 

For  a  discussion  of  the  effect  of  choice  of  the  anchor  equation  see  [30] . 

Hopf  Bifurcation 

The  Hopf  bifurcation  points  are  located  by  finding  a  solution  < v»  >  on  the  steady 
state  branch  where  the  eigenvalues  of  fu(v,XQ)  are  Purely  imaginary.  The  periodic 
solution  separating  itself  from  the  steady  state  branch  at  these  points  is  describable  by 
j  the  following  asymptotic  expansions  [37,  38] 


u(t,e)  -  v  ♦  e»(t)  +  0(e  ) 


P ( £ )  -  PQ  +  0(e  ) 


*0  ♦  0(e) 


(15) 

(16) 
(17) 


where  ±w/p  are  the  eigenvalues  at  the  Hopf  point  and  where  e  is  the  amplitude  of  the 


orbit  ♦(t)  is  the  corresponding  2*-periodic  solution  of  the  linearized  problem 


(18) 


(15)  through  (18)  suggest  the  following  procedure  for  switching  onto  the  periodic  solution 
branch.  Assume  the  steady  state  branch  has  been  computed  so  that  the  Hopf  points  have  been 
located  and  PQ  and  $(t>  computed.  Normalize  $(t)  such  that  * ♦ ■ 2  -  1.  Next  solve 
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rV»3h,  ’•  • 


J 


I 


(12)  but  replace  (14)  with  the  nearly  Identical  equation 


r  T 

A(u)  -  /  (u(t)  -  v)  ♦*(t)dt  -  0 


Under  reasonable  assumptions  the  modified  set  of  equations  can  be  shown  to  be  well-posed 
near  the  Hopf  bifurcation  point,  so  that  computational  procedures  can  be  based  upon  them. 
Discretisation 

The  numerical  solution  of  (12)  requires  the  discretization  of  (10)  and  evaluation  of 
the  integrals  in  (13)  and  (14).  There  are,  of  course,  many  ways  to  do  this  and  to 
illustrate  how  a  typical  computation  might  proceed,  we  describe  in  some  detail  a  specific 
choice  of  the  discretization.  The  differential  equation  is  approximated  by  the  method  of 
collocation  at  Gauss  points  with  piecewise  polynomials  that  belong  to  the  class  C(0,2*) 
[39-41].  In  our  program  the  number  of  collocation  points  per  mesh  Interval,  and  hence  the 
method  order,  can  be  selected.  For  simplicity  we  describe  the  technique  for  the  specific 
choice  of  two  collocation  points  per  mesh  interval.  Introduce  a  mesh 
(0  -  tQ  <  t1  <  ...  <  tN  -  2*},  At.  5  t^  -  t^,  (1  <  j  <  N) ;  and  for  each  j, 

introduce  the  Lagrange  basis  polynomials 

wj_1(t)  -  2(t  -  tj_1/2)(t  -  t^/At*  , 

WJ,-V2<t)  “  "4(t  ’  tj-1)<t  ‘  tj>/Atj  ' 


WJ.0!t>  ‘  2(t  ’  VlUt  "  Vl/2>/Atj  ' 


tj-1/2  “  2  ltJ-1  *  tJ)  ' 


Relative  to  each  subinterval  [t  ,t  ],  the  Gauss  points  are  C  »  t.  ,  -  At  /(2/J) 

J  1  J  j/'  J“V2  j 

end  2  *  tj-1/2  +  The  collocation  method  now  consists  of  finding 

Pj(t)  ‘  Wj,-1(t)Uj-1  *  W3.-1/2(t)Vl/2  *  Wj.0<t,Uj'  1  4  1  *  *  ' 


such  that 


i  -1,2»  1  <  j  <  N 


(20) 


and  (0 )  -  Pn(2»),  i.e.. 


With  the  above  basis  selection  u^ 
the  continuous  problem  at  t  ^  and 
discretized  by  Simpson's  rule  giving 


U0  ‘  “n  • 

and  are  to  aPProx^lute  the  solution  u(t) 

*"j-1/2  respectively.  The  integral  in  (13)  can  be 


(21) 


of 


91<x  -  V2  +  92<P  -  V2  + 


(22) 


Vj-1>2  ♦  «Vl/2  '  Vl/2)2  +  <Uj  _  V2"6 


2 

a  , 


(Accuracy  is  not  actually  important  in  this  particular  equation,  since  it  does  not  affect 
the  overall  order  of  the  numerical  scheme.)  The  anchor  condition  (14)  is  also  discretized 
by  Simpson's  rule  to  give 


N  At 


^  t<uj-i  ’  vl-i)Tf(VV  + 


(23) 


+  4(Uj-1/2  -  Vj-1/2)Tf(Vl/2'V  *  <Uj  -  Vj>Tf<vJ'X0))  ■  ° 


except  at  a  Hopf  bifurcation  point,  where  we  replace  (23)  by  the  corresponding 
discretization  of  (19).  Equations  (20)-(23)  are  the  discretized  form  of  (12). 

Application  of  Newton's  method  to  equation  (12)  requires  the  solution  of  linear 
systems  of  the  form 
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U  V  V 

N  N  N, 
u  P  A 


A  0  0 

u 


Ap 


- 


,  V  -  0,1,2,. 


(24) 


Here  C  is  a  matrix  of  dimension  2Nn  by  (2N  +  1)n,  having  block-diagonal  structure, 
u 

Each  block  can  be  subdivided  into  6  square  matrices  of  dimension  n,  that  is, 

v 

w'  , (5 .  .  )I  -  zrz  w,  .(£.  )f  (p^lS.  .  )»A  )  , 
j,i  ],k  n  2*  3,i  3,k  u  j,k 

where  I  is  the  n  by  n  identity  matrix,  and 

Pj(t)  “  Wj-1<t)Uj-1  +  W j,-1/2(t)U j-1/2  +  Wj,0(t>Uj  • 

C  and  C,  are  column  vectors  of  dimension  2Nnj  N  and  A  are  row  vectors  of 
p  A  u  u 

dimension  (2N  +  1)nj  while  Np  and  N^  are  scalars. 

With  the  Jacobian  matrix  partitioned,  the  linear  system  to  be  solved  for  each  Newton 
iteration  has  the  form 


A  B 


C  D 


An  efficient  solution  proceeds  as  follows:  First,  decompose  the  matrix  as 

C :)(:  :)■ 

i.e.,  decompose  A  into  A  *  LU,  solve 

LP  -  B 


T  T  T 
0  Q  -  C 
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and  set 


T  -  D  -  QP  . 

Then  find  intermediate  quantities  ^  and  by  solving 


and  letting 


Finally  obtain  the  solution  from 


L*1  -  h1 


*2  =  h2  -  * 


TX2  -  *2 


Ux1  -  -  PX2  . 


Pivoting  is  generally  necessary  to  solve  the  next  to  the  last  equation  with  matrix  T. 

Note  that  we  can  compute  the  LU- decomposition  of  X  very  rapidly  due  to  its  simple 
structure.  The  indicated  algorithm  for  solving  the  linearized  system  is  not  necessarily 
the  most  stable  although  it  has  been  used  successfully  in  a  variety  of  problems.  Other 
schemes  should  be  considered  also. 

In  order  for  Newton  iteration  to  converge  to  the  intended  solution  of  (24)  an  accurate 
initial  approximation  must  be  provided.  In  general  this  can  be  done  by  extrapolating  from 
the  two  preceding  points  on  the  solution  branch.  Xt  a  Hopf  bifurcation  point  the 
asymptotic  expansions  ( 15)  —  ( 17>  can  provide  the  initial  approximation. 


RESULTS 

We  now  present  some  of  the  more  interesting  results  obtained  by  using  the  above 
numerical  methods  to  study  the  dynamics  of  a  CSTR  with  A  ♦  B  ♦  C  kinetics  which  arise 
from  varying  the  Damkohler  number.  These  results  exhibit  dynamic  phenomena  which  are  much 
more  complex  than  the  A  *  B  case  and  illustrate  the  utility  of  our  techniques  in  tracing 
entire  periodic  solution  branches  containing  unstable  oscillations  and  multiple  turning 
points. 

For  a  particular  combination  of  system  parameters  [36),  we  summarize  the  multiple 
steady  and  oscillatory  states  on  a  response  curve  with  the  Damkohler  number  as  the  abscissa 
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and  the  I<2  norm  of  the  steady  state  or  periodic  solution  as  the  ordinate.  The  first 
example  is  illustrated  in  Piure  1  which  corresponds  to  the  parameter  values  B  «  9.0, 
a  =  1.0,  6  -  2.0,  and  o  -  0.01.  For  all  values  of  the  Damkohler  number,  the  steady 
state  is  unique  with  exchanges  in  stability  at  the  four  Hopf  bifurcation  points.  The 
application  of  the  switching  and  continuation  procedures  discussed  in  the  previous  section 
enabled  us  to  trace  one  branch  of  periodic  solutions  joining  the  orbits  bifurcating  from 
the  first  and  second  points  and  another  branch  joining  the  orbits  from  the  third  and  fourth 
points.  With  the  notable  exception  of  two  regions  of  oscillatory  states,  this  example  is 
similar  to  the  A  ♦  B  CSTR.  As  D  is  increased  the  norm  of  the  steady  state  solution 
increases  until  the  first  Hopf  point  (D  •  0.243)  is  crossed  and  the  reactor  begins  to 
oscillate.  The  oscillations  cease  as  the  second  point  (D  «  0.29)  is  reached.  Then  the 
steady  state  value  of  lul  increases  with  D  and  the  same  dynamic  pattern  is  repeated. 

The  results  presented  in  Figure  2  are  certainly  not  comparable  to  those  found  in 
simpler  reactors.  In  this  case  we  again  find  a  unique  steady  state  branch  containing  four 
Hopf  points,  but  now  the  orbits  bifurcating  from  the  second  and  third  points  are  joined 
together  as  are  the  solutions  from  the  first  and  fourth  points.  Close  examination  of  the 
inner  branch  reveals  that  both  solutions  bifurcate  to  the  left  of  the  Hopf  point  with  the 
lower  bifurcation  directed  nearly  straight  up  with  a  turning  point  located  close  to  the 
steady  state  branch.  The  periodic  branch  joining  the  first  and  fourth  Hopf  points  is  even 
more  complex.  Again,  both  bifurcations  occur  to  the  left,  and  the  enlargement  in  the 
figure  shows  that  the  branch  actually  contains  three  turning  points  it  D  *  0.1954, 

D  =  0.2003  and  D  =  0.1997.  Our  techniques  were  able  to  compute  all  the  periodic  solutions 
(both  stable  and  unstable)  along  the  branch  and  continue  around  the  turning  points  without 
difficulty.  (Note  the  intersection  of  the  periodic  and  steady  state  branches  at  D  *  0.246 
is  not  a  Hopf  bifurcation  point  but  simply  a  point  where  the  norms  of  the  steady  state  and 
the  periodic  solutions  become  identical)  a  stablility  exchange  does  not  accompany  the 
intersection. ) 

The  location  and  direction  of  the  bifurcations  in  Figure  2  agree  with  the  asymptotic 
results  of  Halbe  and  Poore  [36).  Their  stability  analysis  indicates  that  the  solutions 
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bifurcating  from  the  first  and  third  points  are  unstable  while  the  solutions  splitting  from 
the  other  two  points  are  stable.  Given  these  results  and  assuming  stability  exchanges 
along  the  periodic  branches  occur  only  at  the  turning  points  (preliminary  calculations 
confirm  this  assumption),  we  can  explain  the  complex  dynamic  behavior  exhibited  by  the 
reactor.  Consider  first  the  reactor  operating  at  steady  state  with  D  *  0.23,  and  note  that 
this  state  is  surrounded  by  a  large  amplitude,  stable  limit  cycle  and  a  smaller  amplitude, 
unstable  limit  cycle.  As  D  is  decreased,  we  remain  on  the  steady  state  branch  until  the 
Hopf  point  (D  *  0.2186)  is  crossed  at  which  time,  the  reactor  moves  out  along  the  periodic 
branch  and  begins  to  oscillate.  The  reactor  operates  in  this  periodic  mode  until  the 
turning  point  at  0*0.2172  is  encountered,  and  if  D  is  decreased  below  this  value,  the 
reactor  must  jump  to  the  stable,  larger  amplitude  oscillations  above.  Transitions  from 
periodic  states  at  turning  points  are  often  seen  in  the  A  ♦  B  CSTR  [7)  but  in  that  case 
the  reactor  must  always  jump  to  a  lower,  stable  steady  state. 

Jump  transitions  between  periodic  states  can  occur  in  a  slightly  different  fashion 
along  branch  joining  the  first  and  fourth  Hopf  points.  If  we  further  decrease  D  from 
0.217  to  just  below  0.1997,  the  reactor  will  drop  from  the  upper  solution  to  another  stable 
orbit  of  smaller  amplitude.  The  opposite  effect  is  of  course  achieved  by  now  increasing 
D  past  the  turning  point  at  0.2003  so  that  the  hystersis  phenomena  associated  with 
multiple  steady  states  now  occurs  with  oscillatory  states.  S-shaped  branches  of  periodic 
solutions  are  not  reported  in  earlier  studies  and  could  explain  similar  periodic 
transitions  reported  in  the  tubular  reactor  with  an  A  +  B  reaction  [8] . 

The  next  three  figures  locate  the  concentration  and  temperature  at  various  points 
along  the  solution  branch  and  plot  the  variables  as  time  ranges  from  0  to  2*.  As  the 
reactor  moves  around  the  branch,  the  amplitude  of  the  solutions  becomes  larger  and  the 
temperature  oscillation  becomes  increasing  complex. 

Figure  6  illustrates  a  response  curve  containing  again  four  bifurcation  points  and  now 
a  region  of  three  steady  states.  The  first  two  Hopf  points  are  tied  together  by  a  solution 
branch  like  those  in  Figure  1,  but  the  orbits  emanating  from  the  third  and  fourth  points  are 
not  connected.  In  the  case  of  the  third  Hopf  point,  the  period  of  the  oscillations 
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Figure  5.  Temperature  vs.  Time 


increases  as  one  moves  out  along  the  branch  and  becomes  infinite  as  the  branch 
terminates.  The  next  figure  illustrates  the  increasing  amplitude  of  the  limit  cycles  as 
the  orbits  apparently  approach  a  separatrix  in  a  complete  phase  plane.  (He  have  not  traced 
all  the  branches  to  completion  but  simply  illustrate  the  utility  of  the  numerics  in  such 
computations.  For  example,  the  period  of  solutions  at  the  third  Hopf  point  is  13.81  which 
increases  to  175.64  at  the  last  point  shown  on  this  branch.)  The  occurance  of  'short* 
branches  which  bifurcate  very  near  a  steady  state  limit  point  and  contain  solutions  whose 
period  approach  infinity  are  found  quite  often  in  this  problem.  A  recent  theoretical 
analysis  of  this  particular  type  of  bifurcation  phenomenon  is  presented  by  Keener  [42] . 

The  example  in  Figure  8  shows  a  steady  state  branch  containing  a  region  of  three 
solutions  and  also  three  Hopf  bifurcation  points.  Periodic  reactor  operation  is  possible 
for  a  limited  region  of  D  between  the  second  and  third  Hopf  points  and  for  an  extremely 
narrow  range  of  D  between  the  first  bifurcation  point  and  the  turning  point  contained  on 
the  first  periodic  branch. 

The  results  in  Figures  9  and  10  are  similar  with  the  first  diagram  exhibiting  five 
steady  states  and  two  branches  of  oscillatory  solutions  which  do  not  connect.  (The  steady 
state  branch  which  crosses  the  raised  x-axis  in  Figure  9  continues  down  to  a  limit  point  at 
D  -  0.062  and  lul  ■  1.25  and  then  turns  back  to  the  origin.)  The  response  curve  in  the 
next  figure  again  contains  two  separate  periodic  branches  as  the  steady  state  branch  is 
stretched  into  a  1-3- 1-3-1  pattern.  The  most  interesting  result  in  this  case  is  the 
direction  of  the  periodic  branches  toward  the  steady  limit  points.  The  periodic  branches 
(presumably  stable)  appear  to  become  nearly  tangent  to  the  limit  points  insuring  the 
existence  of  a  stable  solution  for  all  values  of  D.  As  in  all  other  examples  of  this  type 
of  branch,  the  period  of  the  oscillation  becomes  increasing  large  as  we  move  away  from  the 
bifurcation  point.  Analytical  results  concerning  the  above  phenomenon  have  recently  been 
obtained  by  Golubitsky  and  Langford  [43,  44). 

The  last  figure  shows  a  1-3-1  steady  stats  curve  with  two  periodic  branches.  The 
second  branch  in  this  example  bifurcates  to  tha  right  but  then  sharply  turns  back  to  the 
left  at  D  *  0.136  creating  of  course  possible  jumps  between  oscillatory  states  and  high 
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Response  Diagram  Illustrating  1-3-1-3-1  Multiplicity  with  Two  Periodic 
Solutions  (B  =  8.0,  a  =  1.0,  6  =  1.0,  c  =  .01) 


temperature  (or  conversion)  steady  states.  While  the  existence  of  turning  points  on 
solution  branches  is  often  seen  in  the  CSTR  case,  they  have  not  been  reported  earlier  on 


branches  which  do  not  join. 

CONCLUSION 

We  present  a  simple  continuation  technique  for  branches  of  periodic  solutions.  This 
technique  is  easy  to  numerically  implement  and  allows  computation  to  proceed  past  turning 
points  as  well  as  along  unstable  branches.  Without  essential  change  the  procedure  can  be 
used  to  switch  solution  branches  at  the  Hopf  bifurcation  points.  The  techniques  were 
applied  to  the  CSTR  with  consecutive  A  ♦  B  ♦  C  reactions.  The  computational  results 
reveal  several  types  of  dynamic  behavior  not  reported  in  earlier  studies  of  this  reactor  or 
the  simpler  A  ♦  B  case. 

It  is  to  be  emphasized  that  the  numerical  methods  presented  here  are  completely 
general  for  nonlinear  autonomous  systems  of  differential  equations.  Extension  of  the 
method  to  more  general  operator  equations  than  those  discussed  in  this  paper  is  not 
difficult,  in  principle.  We  believe  these  techniques  can  find  wide  application  in  the 
bifurcation  analysis  of  mathematical  model'  in  several  fields  of  research. 
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